A MATRIX FRAMEWORK FOR THE SOLUTION OF ODEs: 
INITIAL-, BOUNDARY-, AND INNER- VALUE PROBLEMS 

MATTHEW HARKERt AND PAUL O'LEARYt 

Abstract. A matrix framework is presented for the solution of ODEs, including initial-, bound- 
ary and inner-value problems. The framework enables the solution of the ODEs for arbitrary nodes. 
CO , There are four key issues involved in the formulation of the framework: the use of a Lanczos pro- 

cess with complete reorthogonalization for the synthesis of discrete orthonormal polynomials (DOP) 
C " 3 i orthogonal over arbitrary nodes within the unit circle on the complex plane; a consistent definition 

^\1 ' of a local differentiating matrix which implements a uniform degree of approximation over the com- 

plete support — this is particularly important for initial and boundary value problems; a method 
of computing a set of constraints as a constraining matrix and a method to generate orthonormal 
^^ , admissible functions from the constraints and a DOP matrix; the formulation of the solution to the 

'^l^ 1 ODEs as a least squares problem. The computation of the solution is a direct matrix method. The 

worst case maximum number of computations required to obtain the solution is known a-priori. This 
makes the method, by definition, suitable for real-time applications. 

The functionality of the framework is demonstrated using a selection of initial value problems, 
Sturm-Liouville problems and a classical Engineering boundary value problem. The framework is, 
however, generally formulated and is applicable to countless differential equation problems. 
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1. Introduction. There are a number of papers in which the Taylor Matrix 
is used to compute solutions to differential equations [TTl [Tl]. These methods use 
►^ ' the known analytical relationship between the coefficients s of a Taylor polynomial 

^SJ ' E^nd those of its derivatives s to compute a differentiating matrix D. The matrix D 

together with the matrix of basis functions arranged as the columns of the matrix B 
en , are used to compute numerical solutions to the differential equations. The method of 

the Taylor matrix was also extended to the computation of fractional derivatives [12] . 



■"nI" i The problem associated with this approach is that the computation of the numerical 



solutions requires the inversion of the Vandcrmonde matrix, a process which is known 
to be numerically unstable, and dependent on the degree and node placement. The 
advantage of the Taylor approach lies in its ability to yield a solution for arbitrary 
nodes. 

^^ I A Chebyshev matrix approach was presented by Sezer |26j . The approach is fun- 

^ . damentally the same as for the Taylor matrix, whereby the Chebyshev polynomials 

are used as an alternative to the geometric polynomials. The main restriction associ- 
ated with the Chebvshev polynomial approach is that the numerical solution to the 
differential equations is restricted to the locations of the Chebyshev points; this lacks 
the generality needed for many differential equations and applications. 

Podlubny introduced a matrix approach to discrete fractional calculus \20\ and 
later extended this work to partial fractional differential calculus [221 [21] . Triangular 
strip matrices play a central role in the work; they are used to perform integration. 
They implement the integration from a lower to an upper bound (or vice versa), 
whereby the errors accumulate as the integration proceeds. This poses a problem if 
inverse problems are addressed, since it gives the solution an implicit direction and a 
different accumulation of errors if the problem is solved from lower to upper bound 
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or from upper to lower. Furthermore, it is assumed that the initial value is zero. This 
makes the method unsuitable for arbitrary boundary conditions. An early source of 
this formulation was proposed by Courant et al. [5 , (a later English translation of 
the paper is available ^6j). 

A matrix solution specific to Sturm-Liouville problems was presented by Amodio [2] . 
The method is specifically restricted to Sturm-Liouville problems; furthermore, it only 
supports solutions on regularly spaced nodes. The results are correspondingly modest 
for problems where the Chebyshev points yield better solutions, e.g., in the solution 
of the truncated hydrogen equation. A number of matrix approaches based on the 
Numerov method, and modifications of this method, have also been presented [15] for 
the solution of Sturm-Liouville problems, however, these methods can not be extended 
to ODEs in general. 

In this paper we formulate a general matrix framework for the solution of ordinary 
differential equations, with arbitrary initial-, boundary-, or inner values, the main 
contributions of the paper are: 

1. The proposal of a consistent framework of matrices and solution approaches 
which can be applied to initial-, boundary-, and inner-value problems; 

2. The implementation of new approached to the synthesis of discrete orthonor- 
mal basis functions, with and without weighting; 

3. Generating differentiating matrices which are of constant degree of approx- 
imation over the complete support. It is particularly important that the 
degree of approximation is consistent at the ends of the support if initial and 
boundary value problems are to be solved satisfactorily; 

4. The derivation of a means of synthesizing constrained basis functions which 
form orthonormal matrices. This basis functions span the space of all solu- 
tions which fulfil the constraints. They can be used as admissible functions 
in a discrete equivalent of a Rayleigh-Ritz method; 

5. The formulation of the solution of the ODEs as least squares approximations. 
In this manner there is no accumulation of errors. 

This paper is structured as follows: In Section ^ the framework for the gener- 
ation of all the matrices required to formulate differential equations as matrix linear 
differential operators is presented. Section ([3]) presents the approach to discretization 
of the differential equations and their solution as a least squares minimization is pre- 
sented. The required conditions for a unique solution are derived and two solution 
approaches are presented: a direct solution in the case of a unique solution and the 
implementation of a discrete Rayleigh-Ritz method for eigenvalue/eigenvector solu- 
tions, e.g., as encountered in the solution of Sturm-Liouville problems.. Finally, in 
Section Q the performance of the proposed framework is tested with a series of initial- 
value problems, Sturm-Liouville problems and a classical Engineering boundary value 
problem. 

2. Algebraic Framework. In this section we derive the structure and methods 
for the synthesis of all matrices required for the discretization and solution of ordinary 
differential equations. 

2.1. Quality Measure for Basis Functions. An objective measure for the 
quality of a set of basis functions is required if the sources of numerical error are 
to be determined and the best synthesis method is to be selected. In this paper 
continuous polynomials are considered which form orthogonal bases when evaluated 
over a discrete measure. The basis functions 6^, i.e., the polynomials evaluated at 
discrete points, can be concatenated to form a matrix, B = [bi . . .6„]. The discrete 



orthogonal polynomials (DOP) are characterized by the relationship, 

B'^WB = I, (2.1) 

where W is the weighting matrix. The Gram matrix is defined as G = B""" W B. 
Consequently, the orthogonal complement G-*" = I — B""" W B = should be a matrix 
containing only zeros. However, this is not the case, due to the loss of orthogonality 
in the three term relationship resulting from numerical errors. These numerical errors 
determine the quality of the basis functions and for which we require a measure. The 
determinant of G has in the past been used as a measure for the quality Eg — det G 
of the basis functions. However, this measure does not yield stable estimates [H 
Chapter 2, Sec. 2.7.3]. We propose the Frobenius norm of G-*^ as an error measure, 
i.e., €f = IIG^IIi;-, this is the sum of the square of all errors w.r.t. the orthogonality 
of the basis functions, ep > 0. This is a posteriori measure, i.e., we compute the 
basis functions and then determine their quality. Wilkinson 1271 points out that a- 
priory prediction of error bounds yield unreliable results and a posteriori analysis is 
preferred. The numerical results obtained for different synthesis procedures can be 
found in Section (14.11). 



2.2. Numerically Stable Synthesis of Basis Functions and their Deriva- 
tives. Gram [5] proposed what is now known as the Gram-Schmidt orthogonalization 
process to generate polynomials 4 . The Gram-Schmidt process is, however, numer- 
ically unstable [8] Chapter 5] and errors accumulate as the number of integrations 
increases, i.e., with increasing polynomial degree. This precludes the synthesis of 
polynomials of higher degree with this method. Considerable research has been per- 
formed on discrete polynomials and their synthesis [16l[29l|30l[28l[l0l[3Tll32l|3]. The 
research was primarily in conjunction with the computation of moments for image 
processing. None of these papers present a method which is capable of synthesizing 
discrete orthogonal polynomials of high quality for arbitrary nodes located within the 
unit circle on the complex plane. 

Here it is proposed to synthesize the polynomial basis functions using a Lanczos 
process with complete reorthogonalization |8, Chapter 9, p. 482],J/7j. The procedure 
can be summarized as follows: Given a vector x of n nodes with mean x, i.e., the 
points at which the differential equation is to be solved: first compute the two basis 
functions bo, bi and initialize the matrix of basis functions B, 

bo = l/V^. bi= ,1 '^~!',| and B = [bo,6i]. (2.2) 

The remaining polynomials are synthesized by repeatedly performing the follow- 
ing computations: 

1. Compute the polynomial of the next higher degreqj, 

6„=6io6„_i; (2.3) 

2. perform a complete reorthogonalization, 

b„ = 5„ - B B^ b„ (2.4) 

= {I-BBT}6„ (2.5) 



^The symbol o represents the Hadamard product. 



by projection onto the orthogonal complement of all previously synthesized 
polynomials. It is important to note that the reorthogonalization is w.r.t. to 
the complete set of basis functions, not just the previous polynomial. 

3. Normalize the vector, 

bn = jr^, (2.6) 

4. and augment the matrix of basis functions, 

B=[B,b„]. (2.7) 

This procedure yields a set of orthonormal polynomials from a set of arbitrary nodes 
located within the unit circle on the complex plane. Although in [7] the Lanczos 
process is used to compute discrete orthogonal polynomials, the authors seem to have 
overseen the possibility (necessity) of using complete reorthogonalization at each step 
of the polynomial synthesis. 

By taking the derivative of the recurrence relationship w.r.t. x, we obtain the 
equations required to simultaneously synthesize the differentials of the polynomials. 
This procedure appears in |13| for the Legendre and Chebyshev polynomials. Here 
the method is generalized to the synthesis of polynomials from arbitrary nodes. With 
this, the synthesis procedure delivers a set of orthonormal basis functions B and their 
derivatives B. 

2.3. Weighted Basis Functions. A set of discrete basis functions in matrix 
form, Bu, are orthogonal with respect to a weighting matrix W if, 

b;i;wb^ = I. (2.8) 

In the case of a weighting function w{x) the weighting matrix is given by W = 
diag {w {x i) .. .w{xn)}- Given a set of orthonormal basis functions B and a posi- 
tive definite weighting matrix W, there exists a set of weighted basis functions B^^, 
such that B„ = B U, whereby U is a full rank upper triangular matrix. Substituting 
into Equation (|2.8p yields, 

U^B^WBU^I. (2.9) 

Since U is full rank, we may invert it to obtain, 

B'^WB^ U^'^U"^ (2.10) 

The Cholesky decomposition chol {A} of a matrix exists and is unique such that 
A = G G""" if A is real positive definite. The matrix G is a full rank lower triangular 
matrix. Consequently, the Cholesky decomposition chollB""" WB} exists if W is real 
positive definite, since B is orthonormal. Applying the decomposition yields, 

B^WB = GG'^ = U"'^U"^ (2.11) 

The sought matrix U is clearly given by, 

U = G"^. (2.12) 

With this the weighted basis functions are fully defined. The condition number of the 
basis functions depends soley on the condition number of the weighting matrix W. In 
the case where the weighting matrix is derived from a weighting function w{x)^ the 
condition number is determined by the extreme values of w{x). 



2.4. Differentiating Matrices. There are both global [TTl [H HllSSl IS] and 
local [SJ [15J ISl [ISl [22] approaches to computing discrete estimates for derivatives. 
Global methods proposed in the past have used the known relationship between the 
coefficients of a polynomial and the coefficients of the derivative of the polynomial to 
compute a differentiating matrix. 

The computation of a differentiating matrix from polynomial bases proceeds as 
follows: The spectrum of the signal y with respect to the basis functions B is computed 
as, 



B+y. 



(2.13) 



For example, B may be the Vandcrmonde matrix; this is case with Taylor methods |lli 
[T¥[ FT^ . The relationship between the spectrum s and the spectrum of the derivatives 
is given by. 
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(2.14) 



Consequently, 



y = BW\B+ y=Dy. 



(2.15) 



That is, the differentiating matrix is computed as, D = BMB+. In the case of 
the Taylor (Vandermonde) matrix this involves computing the pseudo-inverse of the 
Vandcrmonde matrix: with all the associated numerical problems. In the case of the 
Chebyshev polynomials [221 [13] B"*" — B^ and a different matrix M is required, see [2S] 
for details. The method is not appropriate if arbitrary nodes are required, e.g. this 
may be required if the framework is to be used to solve the problems associated with 
monitoring mechanical structures [19] . The advantage of Global methods is that they 
deliver a differentiating matrix which is valid for the complete support. 

The solution chosen here is to compute D form the basis functions and their 
derivatives, i.e., given B and B an appropriate derivative operator, D, should have the 
property that. 



Post-multiplying by B""" yields 



DB = B. 



Dr = DBB'^:= BB'^. 



(2.16) 



(2.17) 



If the basis function set is complete, i.e., B B"^ then the above equation yields the 
differentiating matrix directly, 



d = bbJ 



(2.18) 



This computation is valid for arbitrary nodes. If a truncated, i.e., an incomplete, set 
of basis functions is used then B B""" is the projection onto the the basis functions B 
and Db is then a regularizing differentiating matrix. The matrix Ds can be applied 
to the computation of estimates for derivatives in the presence of noise. 
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Fig. 2.1. Rank deficiency of the differentiating matrix D as a function of the degree of a Gram 
polynomial, when D is synthesized using Equation i2.18[ l. 



A differentiating matrix should be rank-1 deficient; it should have the constant 
vector as its null space, i.e., Dla = 0. The properties of the D are, however, depen- 
dent on the nodes being used, e.g. the Chebyshev nodes permit global differentiating 
matrices for very high degrees. With other sets of nodes the condition number of a 
differentiating matrix can increase with the degree of the polynomial being used. At 
some point the matrix starts to have additional null spaces, which are associated with 
numerical errors occuring due to insufficient numerical precision, this effect is shown 
in Figure (|2.ip for the Gram polynomials. 

Once the condition number of a global differentiating matrix has degenerated 
below an acceptable level, it becomes necessary to compute local approximations [251 
[T8] . Courant [5l|6] proposed, in 1927, using both forward and backward differences to 
compute estimates for the first derivative. The method has also been used in [501 112] 
to this end. More commonly the tri-diagonal matrix, shown here for 6 points, 
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is used to compute a discrete local estimate for the first differentiao. This operator 
is only of degree d — 2 accurate in the core of the approximation and at both ends 
of the support only of degree d = I accurate. This makes this discrete operator 
unsuitable for the computation of derivatives at the end of the support, as is required 
for BVPs and IVPs. It is also not suitable for systems whose solutions are locally of 
degree higher than d> 2. Furthermore, this operators assumes equally spaced nodes. 
In [2] higher order finite difference schemes are proposed with end-point formulas. 
However, only equidistant spaced nodes are considered. For example, the appropriate 
three point operator for the Gram Dq^^ and for the Ghebyshev Dc,3 nodes, are. 
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^This is the matrix embedded in the Matlab function gradient. 



and 
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both computed for n — 6 points in the interval — 1 < x < 1. Note the three points 
formulas at both ends of the support. 

In keeping with the formulation of the basis functions for arbitrary nodes: the 
method for local differential approximation is also formulated here for arbitrary nodes. 
A generalized formulation of local differentiating matrix requires the vector a; of n 
arbitrarily placed nodes, the support length Is and the degree d of the approximation. 
Only odd support lengths are considered here, to avoid the need for forward and 
backward formulas. It is convenient to define the support length Is = 2ws + 1 in 
terms of the half-width Ws- The vector x of nodes is segmented into m — n ~ 2ws 
overlapping segments, for each segment. 



s{i) = x{i — Ws : i + Ws) Vi G [ws + l,n — Ws 



(2.22) 



a local set of basis functions B^ and derivatives of the basis functions Bs are com- 
puted. Then the differentiating matrix associated with the segment is determined 
Dg = Bs BJ. The first and last segments yield the end-point formulas as required. 
The remaining segment yields the required central formula of coefficients to locally 
approximate the derivative. Clearly, for the inner-segments it is only necessary to 
compute the center row vector of the local differentiating operator D^. The use of 
approximating or interpolating polynomials leads to the generation of differentiating 
matrices with and without regularization respectively. The Wilkinson diagram for the 
general structure of a local differentiating matrix D^ is shown in Equation (j2.23|) for 
the example oi Is = 5 and n = 10. The specific entries in the matrix are a function of 
the spacing of the nodes. 
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(2.23) 



All computations of the local derivative are of length Is and of constant approximation 
degree da — 2ws over the complete support. This is important if derivatives are to be 
computed at the ends of the support; furthermore, errors at the end of the support 
associated with inconsistent approximations will propagate through the entire solution 
when D is being used in the solution of differential equations. This procedure proposed 
here delivers a local differentiating matrix for arbitrary nodes. 



2.5. Defining Constraints. In Section (|2.4I) it was shown that a discrete ap- 
proximation to differentiation can be computed as a linear matrix operator. Con- 
sequently, both differential and integral constraints are linear. In the framework 
proposed here, a constraint is implemented by restricting a linear combination c^ y 
of the solution vector y to have a scalar value d, i.e., 

c^y = d. (2.24) 

This is a very general mechanism, since any constraining function can be implemented 
at a point Xi for which a linear n point expansion around this point exists. To give an 
example, consider the C^ continuous periodicity constraint y{0) = y(l), y(0) = y{l) 
and y{Q) = 2/(1): given the differentiating matrix D and D^, the three constraints can 
be formulated as: 

[l,0,...,0,-l]?/ = cTy = 0, (2.25) 

{D(l,:)-D(end,:)}y = cjy = 0, (2.26) 

{D2(1, :) - D^iend, :)} y = cj y = 0. (2.27) 

Given a set of m constraints, the constraining vectors c^ are concatenated to form 
the matrix C = [ci . . . c„i] and the corresponding scalars di form the vector d = 
[di . . .dm], so that, 

Cy = d. (2.28) 

2.6. Homogeneously Constrained Admissible Functions. Starting from 
a set of basis functions B such that B'^ W B — I , we wish to derive a method of 
synthesizing a set of constrained basis functions Be which fulfil the conditions: 

BJWBc = I, C'^Bc = and B^ = BX, (2.29) 

i.e., the constrained basis functions form an orthonormal basis set with respect to 
the weighting matrix W. If B is orthonormal, i.e., B""" B = I then so is Be. The 
constrained basis functions fulfil the homogeneous constraints defined by C. If B is 
complete then it spans the complete nx n space, given p — rank (C), i.e., the number 
of independent constraints, B^ is of dimension n x (n — p) and spans the complete 
space in which the constraints are fulfilled. Consequently, all possible vectors y which 
fulfil the constraints are given by, 

y=Bca (2.30) 

where a is an n — p vector. 

A solution to the task of determining X was presented in [19]; however, a more 
succinct derivation is provided here. The conditions from Equation (|2.29p require, 

C'^BX = (2.31) 

and with this X must lie in the null space of C''" B. Applying QR decomposition to 
BT C yields, 

QR = B'^C, (2.32) 

and consequently, 

X^QR = (2.33) 



The matrices Q and R are partitioned according to the span and null space of B""" C, 



Q = [Q„Q„] and R 



(2.34) 



with Rs of dimension p x p. The n x p matrix Q^ forms a basis set for the span and 
the n X [n — p) matrix Q„ forms a basis set for the null space of B C. Consequently, 

X^Q, =0 and (X^ Qn)'^ WX^ Q„ = I. (2.35) 

Now applying an RQ decomposition to Q„ yields, 

RQ„ = Q„. (2.36) 

R is orthonormal, since both Q„ and Q„ are by definition orthonormal. Now, selecting 
X = R yields X""" R Q„ = Q„, and with this all the conditions from Equation (j2.29p are 
fulfilled. The matrix X being orthonormal ensures that Be fulfils the same orthonormal 
condition as does B. Furthermore, X has an implicit partitioning. 



X = 



Xi 

X2 



(2.37) 



whereby, Xi is a p x (n — p) block matrix and X2 is a (n — p) x (n — p) upper triangular 
matrix. This structure ensures that the number of roots in the constrained basis 
functions Be is ordered in the same manner as in B. 

3. Discretizing and Solving Ordinary Differential Equations. In the pre- 
vious section all the matrices required for the discretization of ordinary differential 
equations were derived. In this section the discretization of initial-, boundary- and 
inner value problems is presented together with the associated methods of solving the 
resulting matrix equations. 

3.1. Initial Value Problems. In this paper we are considering the solution of 
linear ordinary differential equations with constant or variable coefficients, they can 
in general be formulated as, 

pu{x)y^''\x) . . . + pi{x) y^^\x) + po{x) y{x) = g{x) (3.1) 

to which a set of k constraints are required to ensure a unique solution. The term 
yW(^x) represents the k^^ derivative of y{x). Given the matrices derived previously, 
the discretization of Equation (|3.ip is direct and simple, each term pk{x)y'^^\x) is 
dicreteized as follows: The matrix Pfc is formed such that Pfc = diag {pk{x)}, whereby 
Pk{x) is the vector of values obtained by evaluating the function Pn{x) at the vector of 
points x\ the term y'^^\x) is discretized as D*^ y, i.e., the fc**^ power of D, which is the 
differentiating matrix derived in Section (|2.4p . Summarizing, each term is discretized 
as follows, 

Pk{x)y'^''\x)-^PkD''y. (3.2) 

and the vector g — g{x). Applying this to all terms in Equation p.ip yields, 

PfcD*^^y... + PiDy + Poy = g (3.3) 
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The matrix equivalent ol the linear differential operator L is now defined as, 

L ^ Pfe D^ . . + Pi D + Po, (3.4) 

and the set of k constraints are implemented as defined in Section (12.51) , yielding 

L y = g given C^ y = d. (3.5) 

the matrix C has the dimension n x k. 

A unique solution to the ODE exists only if 



rank 



L 



= n (3.6) 



i.e., the linear differential operator and the constraints must form a full rank system 
of equations. There are many Engineering application where this is not the case, e.g. 
the equations for the vibration of a beam, and Sturm-Liouville problems. A different 
solution approach is proposed for this class of problems in Section p.2[) . 



3.1.1. Solution as a constrained least squares problem. The formulation 
of determining y from Equation (|3.5|) as the solution of a least squares minimization 
problem yields, 

min||Ly-g||2 given C y = d. (3.7) 

y 

This is the well known problem of least squares with equality constraints (LSE). 
Efficient and accurate solutions can be found in [U Chapter 12]. This method will 
yield solutions for ODEs with consistent constraints and a least squares solution in the 
case of over-constrained systems and perturbed systems. It is not a suitable approach 
for Sturm-Liouville type problems. 

The worst case number of floating point operations (FLOPS) required to perform 
the computation is known a-priori. This, by definition, makes the method suitable 
for real time applications. 

3.1.2. Spectral Reqularization. Spectral regularization is introduced here to 
limit the number of zeros in the basis functions and with this to reduce the errors 
associated with aliasing. Assuming y can be sufficiently accurately approximated by 
a series of r orthonormal basis functions, we may write, 

y = B,a, (3.8) 

whereby B^ — B(:, 1 . . . r). Now defining L,. = L B,. and C,. = B^ C, and substituting 
into Equation (|3.7p yields, 

min ||L,. a — g||2 given Cj a = d, (3.9) 

a 

whereby the series coefficients a are to be determined. In addition to introducing 
regularization, the truncated basis functions also reduce the size of the LS problem 
to be solved. 
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3.1.3. Solution of Homogeneously Constrained IVPs. Homogeneously Con- 
strained IVPs for a special subclass of problems for which there is a particularly simple 
solution. Let the solution y be a linear combination of a set of constrained basis func- 
tions, i.e., y = BcQ!, which fulfil the homogeneous constraints C""" Be = associated 
with the IVP. Equation (13.71) now simplifies to the unconstrained least squares prob- 
lem, 



min||LBca-g||2. (3.10) 



The solution of which is. 



a = {LBJ+ff (3.11) 

since null {L Be} = if a unique solution exists. Consequently, 

y = Be{LBJ+g (3.12) 

3.2. Sturm-Liouville and Boundary Value Problems. A Sturm-Liouville 
problem is a second order ODE with the following structure. 



_d_ 
dx 






+ g{x)y = Xw{x)y, (3.13) 



in the finite interval xi < x < a;„, where p{x), g{x) and w{x) are real-valued strictly 
positive. Additionally there are two boundary conditions which are most commonly 
formulated as, 

aiy{xi)+a2v{xi) ^Q, (3-14) 

hiv{x2) + b2y{x2)=Q. (3.15) 

There are some important properties of Sturm-Liouville equations |15j which must be 
considered when implementing a discrete solution: 

1. All eigenvalues are real and there is no largest eigenvalue, i.e., there are an 
infinite number of eigenvalues and A™ — >■ cx) as to — >■ oo. Given a set of n 
discrete points x there can theoretically only be n eigenvalues; 

2. The m}^ eigenf unction has ra zeros on the interval a < a; < 6. However, 
given n points (samples) only functions with a maximum of n/2 zeros can be 
discribed without aliasing. Consider the Sturm-Liouville equation y — Xy = 
with the constraints y{0) — and y{Tr) — 0. This equation is known to have 
the eigenfunctions ^„i (x) = v2 sin (to, tt x) . Consequently, a discrete solution 
can only model the first n/2 eigenpairs correctly. 

3. The eigenfunctions are orthogonal with respect to the weighting function 
w{x), i.e., /^ w{x)<i>i{x)<i>j{x) = S{i,j). ^^ 

The general Sturm-Liouville problem formulated in Equation (I3.13|) with its cor- 
responding boundary conditions can be discritized directly as, 

{DPD-G} y = -AWy given C"^ y = 0. (3.16) 

whereby, P = diag{p(a;)}, G = dia.g{g{x)} and W — dia.g{w{x)}. A direct solution 
of this equation will, however, yield unstable results due to aliasing. 

We now introduce a set of weighted and constrained basis functions B^, which fulfil 
the orthogonality condition B^ W Buj = I and boundary conditions C"^ B = 0. These 
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basis functions are admissible functions for the Sturm-Liouville problem. The number 
of zeros in the basis functions increases from left to right in the matrix. The number 
of zeros in the admissible functions is limited, so as to avoid aliasing, by truncating 
to the first k = n/2 basis functions, i.e., B^ = Bi„(:,l : k). The eigenfunctions 
are now found as linear combinations of these admissible functions, i.e., y = B^ ct. 
Substituting this into Equation p.l6|) yields, 

{DPD-G} BaQ; = -AWBaa. (3.17) 

Pre-multiplying both sides by Bj now yields, 

Bj {DPD-G} B„a = -Aa. (3.18) 

since Bj\N Ba = I. Now defining La — Bj {D P D + G} Ba yields a standard eigen- 
vector problem, 

{la +X\} a = 0, (3.19) 

y=Baa. (3.20) 

Solving Equation (j3.19p for the eigenvalues A^ and the eigenvectors cti, then back 
substituting a^ into Equation p.20p yields the desired eigenfunctions. 

It is important and interesting to note the the matrix L^ is of dimension n/2 x 
n/2, in contrast to the original matrix L = D P D — G which is of dimension n x n. 
Consequently, dealing with the aliasing has also reduced the size of the eigenvalue 
problem to be solved. In the worst case an eigen-decomposition is of complexitj|j 
between 0(n^) and 0{n^). The improvement in speed is then in the range of a factor 
of 4 to 8, while simultaneously improving the accuracy of the solution. However, 
some of the computation gains are spent on additional pre- and post-calculations. 
A consequence of Equation (j3.20|l is that the matrix of eigenvectors a, contains the 
spectrum of the eigenfunctions with respect to the basis functions used, i.e., the 
Rayleigh-Ritz coefficients. 

4. Performance Testing. In this section a selection of examples are presented 
to demonstrate the functionality of the proposed methodQ. 

4.1. Quality of Basis Functions. The first test addresses the quality of basis 
functions, since these form the basis for all subsequent calculations. The following 
polynomials are compared: a set of Gram polynomials generated using Gram-Schmidt 
orthogonalization [^ ; a set of Chebyshev polynomials generated using the recurrence 
relationship |23j ; a Vandermonde matrix and a set of polynomials synthesized using 
the method proposed in this paper. The Frobenius norm of the projection onto the 
orthogonal complement of the Gram matrix is used as an estimate of the total error. 
The number of significant digits is then estimated to be d = — logiQ{ep). Two com- 



putations were performed: Figure (4.1(a) I shows the result for complete polynomial 



sets, i.e., the degree d = n — 1 where n is the number of nodes; Figure (4.1(b) ) is for 
a fixed number of nodes n — 1000 and the degree of the polynomial is progressively 
increased. The results shown in Figure (|4.ip indicate that the algorithm presented in 



•^Indeed there are more efficient algorithms; however, their complexity depends on the structure 
of the matrix and the distance between the eigenvalues. Consequently, no general statements can be 
made about these methods. 

*A MATLAB toolbox D OPbox is available at 

http://www.mathworks.de/matlabcentral/fileexchange/41250 which implements all the methods 
proposed in this paper; furthermore, the code for all the following examples in also provided in 
the toolbox; ODEbox http://www.mathworks.de/matlabcentral/fileexchange/41354l 
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Cliebysliev 



Gram 



Vandermonde 



(a) Complete polynomial basis sets, i.e., d = n- 
1 where n is the number of nodes. 



Chebytihev 




200 300 

Degree 



(b) Number of nodes n = 1000, the degree of 
the polynomial is increasing. 



Fig. 4.1. Estimate of the number of significant digits for different polynomials as a function of 
degree. DOP refers to the discrete orthonormal polynomial synthesis as proposed in this paper. 



this paper generates the most stable sets of polynomials. For this reason the algorithm 
is used for the generation of all bases required in this paper. 

4.2. Initial Value Problems. In each of the following examples the analytical 
solution is compared with the numerical results commuted using the newly proposed 
method and a Runga-Kutta solution with variable step sizqf|. The ode 4 5 is used 
for comparison in all the following initial value problems. 

4.2.1. IVP Example 1. The first example is a second order initial value prob- 
lem with constant coefficients, the equation is [1], 



j/-6y-92/ = 0, with, 2/(0) = 10 and y(0) = -75. 
in the interval < x < 3 The analytical solution to this equations is. 



(4.1) 



(4.2) 



The analytical solution and the results of the numerical solutions are shown in Fig- 



ure (4.2(a)). Non-uniformly spaced nodes were used, with a higher density of nodes 



where y has a higher first derivative. This demonstrates the possibility of generating 
basis functions from arbitrary nodes. The residual errors, see Figure (4.2(b)), with 
the new method is 7 orders of magnitude smaller than with the ODE4 5 method. 

4.2.2. IVP Example 2. The second example is a second order differential equa- 
tion with variable coefficients, the equation is [1], 



2x^ y — xy — 2y = 0, with, y[l) — 5 and y{l) — 0. 
in the interval 1 < x < 10. The analytical solution to this equations is. 



y{x) =x +-^. 
Jx 



(4.3) 



(4.4) 



The comparison of the analytical solution with the numerical solutions using the new 
method and a Runge-Kutta procedure is shown in Figure ()4.3p . The residual error 
with the new method is 3 orders of magnitude smaller that with the Runga-Kutta 
procedure. This example has demonstrated the ability of the proposed method to 
solve initial value problems with variable coefficients. 



'See the MATLAB documentation for details of the ODE45 solver used for this computation. 



14 



Analytical 

o New 

V Runge-Kutta 




10 

a 5 



-5 



(a) Comparison of the analytical solution, the 
numerical solutions using the new method (with 
support length Is = 13, the nodes are placed 
X = 5z^ and z has n = 85 evenly spaced points 
in the interval < z < 1) and using a Runge- 
Kutta procedure. The nodes are shown below 
the curve. 




(b) Residual errors; (top) between the analyti- 
cal solution and the new numerical procedure, 
(bottom) between the analytical solution and 
the numerical solution using the Runga-Kutta 
procedure. 



Fig. 4.2. Computation results for the IVP Example 1 given in Equation 




X 10' 




(a) Comparison of the analytical result, the nu- 
merical solutions using the new method (with 
support length l^ = 13, there are n = 73 evenly 
spaced nodes in the interval 1 < x < 10) and 
using a Runge-Kutta procedure. 



(b) Residual errors; (top) between the analyti- 
cal solution and the new numerical procedure, 
(bottom) between the analytical solution and 
the numerical solution using the Runga-Kutta 
procedure. 



Fig. 4.3. Computation results for the IVP example 2 given in Equation {4-3\ 



4.2.3. IVP Example 3. The third example PQ is a third order non-homogeneous 
differential equation with constant coefficients, the equation is, 

y + 3 y + 3y + y =: 30 e""^ with, y(0) = 3, y(0) = -3, and y{0) = -47 

(4.5) 
in the interval < x <8. The analytical solution to this equation is, 



y{x) = (3-25a;2 + 5x3)e" 



(4.6) 



The comparison of the analytical solution with the numerical solutions using the 
new method and a Runge-Kutta procedure is shown in Figure (|4.4p . Once again 
the residual error with the new method is orders of magnitude smaller that with the 
Runga-Kutta procedure. 
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V Runge-Kutta 
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(a) Comparison of the analytical solution with, 
the numerical solutions using the new method 
(with support length Is = 13, there are n = 73 
evenly spaced nodes in the interval < a:: < 8) 
and a Runge-Kutta procedure. 




(b) Residual errors; (top) between the analyti- 
cal solution and the new numerical procedure, 
(bottom) between the analytical solution and 
the numerical solution using the Runga-Kutta 
procedure. 



Fig. 4.4. Computation results for the IVP Example 2 given in Equation J4-5|). 



4.3. Sturm-Liouville Problems. The test package for Sturm-Liouville Solvers ^ 
has been used as a source of test cases in this sectior^. 

4.3.1. Sturm-Liouville Example 1. The simplest Sturm-Liouville problem is 
chosen as the first example, since the analytical solution is known. This enables the 
investigation of the stability of the numerical computation, i.e., how many eigenvalues 
can be computed to a given accuracy. It is the equation of a vibrating string, 



y + Xy = 0, given y{0) — 0, and y(7r) = 



(4.7) 



in the interval < x < tt. The analytical solution yields the analytical eigenvalues Xk 
and eigenfunctions yk, 



Xk — k, yk ~ sin kx for fc = 1 , 



(4.8) 



The discrete solution has been computed on the interval < a; < tt sampled at the 
corresponding Chebyshev points; however, scaled so that the first and last points lie 
exactly at and tt respectively. For a given set of n points m = n/2 constrained 
basis functions are computed which fulfil the boundary values. These are the admis- 
sible functions used in what is essentially a discrete equivalent of the Rayleigh-Ritz 
method. Two different computations ni = 100 and n2 = 1000, have been performed 
to investigate the behavior of the solution with respect to the number of points used. 
A support length Is = 13 was used during the generation of the differentiating ma- 



trix. The results can be seen in Figure (4.5(a) I and (4.5(b) I respectively. The method 
returns the coefficients of the series of admissible functions required to generate each 
eigenfunction. The matrix of these coefficients is denoted by R. We use the matrix 
S — logj^Q {abs(R)} as a visual representation for the spectrum in the figures, since at 
S(i, j) ~ —16 the numerical resolution of computation environment is reached. This 
makes a visual recognition of when the spectrum fails simple. 

With Til ~ 100 approximately ki — 28 and for 712 = 1000 approximately ^2 = 280 
eigenvalues could be computed with a relative error smaller than 0.1%. This result 



''At this point we feel it is important to note that the framework presented here is generally 
applicable to the solution of ODEs in general and is not a dedicated Sturm-Liuville solver. The use 
of Sturm-Liouville problems as a test cases is to demonstrate this generality. 
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is significantly better than all previously reported results with matrix methods for 
Sturm-liouville problems [TS] . This confirms the numerical stability of the approach. 
It also verifies that the number of eigenfunctions which can be computed to a given 
accuracy scales linearly with the number of nodes used. 




10 20 30 40 

Eigenvector number 71 

(a) Solution for n = 100. 




100 200 300 400 500 

Eigenvector ntunber n 

(b) Solution for n = 1000. 



Fig. 4.5. Solution of the Sturm- Liouville problem corresponding to the vibrating string for 
points. Top: spectrum of the eigenfunctions with respect to the admissible functions {logio{S). 
Bottom: The relative error in the eigenvalue Aj. in %. 



4.3.2. Sturm-Liouville Example 2. The second example is a Mathieu differ- 
ential equation; we have taken this example from i24; Problem 2, with r — 25]. This 
equation arises in the vibration of elliptical membranes. We have chosen this prob- 
lem because it is known to produce a pair of closely located eigenvalues; this should 
enable the test of the resolution of the eigenvalues computed using the proposed 
method. Secondly, the solution to the equation has no known analytical form, this 
makes numerical solutions particularly valuable. The Mathieu differential equation 
is. 



y + 2r cos{2x)y — Xy, with y(0) = 0, and y[n) = 



(4.9) 



in the interval < x < tt. A Chebyshev distribution of n = 1000 nodes are used 
covering the complete interval, a support length Ig = 13 and m = 500 basis functions 
were used for the computation. The result of the numerical computation can be seen 
in Figures (4.6(a) I and (4.6(b)). The pair of expected eigenvalues are computed as, 
Ai = -2.131489^ + 01 and A2 = -2.131486£; + 01. 

The results of these computations are slightly more difficult to interpret abso- 
lutely since the analytical results are not given. It can be said that the expected 
pair of eigenvalues have been found. Secondly, from the nature of the Rayleigh-Ritz 



spectrum and the corresponding computed eigenvalues, see Figure (4.6(b)), suggest 
that approximately the first k — 350 eigenvalue are correctly computed. 

4.3.3. Sturm-Liouville Example 3. This is the truncated hydrogen equation, 
taken from PU Problem 4] . This is an example of s singular Sturm-Liouville equation 
with a limit point non-oscillatory (LPN) end point. Although only one constraint is 
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100 200 300 400 500 
Eigenvector number n 

(a) The first two eigenfunctions of tfie Mathieu (b) Top: Rayleigh-Ritz Spectrum of the Math- 
differential equation, for r = —25; note the neg- icu differential equation. Bottom: the eigenval- 
ative sign for the value of r. ues. 

Fig. 4.6. Solution of the Mathieu differential equation. It is important to note that the coeffi- 
cient r is negative. 



available the equation is well conditioned, since g{x) — >■ (X) as x — >■ 0. Highly accurate 
results for some eigenvalue are available for this equation [24]. These values are used 
to evaluate the accuracy of the computation method presented here. The differential 
equation is, 



y 



1 



y = Xy, with y{0) = LPN, and y(lOOO) = (4.10) 



in the interval < x < 1000. 

The computation was performed using n — 1000 Chebyshev points on the range 
< X < 1000; further, m = 500 basis functions were used, whereby only one constraint 
is applied, i.e., at x = 1000 and a support length of Is = 13. The result of the 
computation are presented in Table (|4.1|) . The known eigenvalues, i.e., Ao, Ag, A17 
and A18 are comparable up the the 10*'\ S***, 6*^ and 5**^ significant digits respectively. 
This indicates a high degree of accuracy, particularly considering that this is a general 
framework for differential equations and not a dedicated Sturm-Liouvalle solver. In [5] 
difficulties with oscillations at the right end point were observed for eigenvalues A^ 
when fc > 8, these difficulties are not observed with the methods proposed here. 





new method 


known value [24 


Rel. Error 


Ao- 


-6.2499999978£;-02 


-6.2500000000£; - 02 


3.4874503285£;-10 


Ag- 


-2.0661156136£;-03 


-2.0661157025£;-03 


4.3009091823£;-08 


Al7 = 


-2.5757218232£;-04 


-2.5757359232£;-04 


5.4741446402£;-06 


A18 — 


2.8740937561£;-05 


2.8739013100^;- 05 


-6.6963370220S-05 



Table 4.1 

Table of eigenvalues for the truncated hydrogen equation, 
obtained using the new method with the known results T2B 



Comparing the numerical results 
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4.4. Boundary Value Problem Example 1. The last test is a classical Engi- 
neering boundary value problem. A cantilever with additional simple support forcing 
the constraint 2/(0.8) = is shown in Figure (I4.7P : this is an example of an inner 
constraint. Furthermore, the system is over constrained, since there are 5 constraints 
placed on a 4*^ order differential equation. It demonstrates the ability of the pro- 
posed framework to solve problems with arbitrarily placed constraints and to solve 
over-constrained systems. The first two admissible and eigenf unctions are shown in 
Figures (4.8(a)) and (4.8(b) I respectively. This computation was performed with 



1/(0) = 



r(i) = o 



1/(0) = 



y(0.8)=0 A i/(l)=0 



Fig. 4.7. Cantilever with an additional support representing an inner-value problem. 
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0.2 0.4 0. 



(a) The first three admissible functions for the 
cantilever with and additional simple support. 



(b) The first three eigenfunctions. 



Fig. 4.8. Admissible functions and eigenfunctions for the cantilever with additional simple 
support shown in Figure j4- 



n — 1001 points evenly distributed along the interval < a; < 1, a total of r = 500 
basis functions were used and a support length Is — 13 was used for the computa- 
tion of the local derivatives. The method presented for this problem is a discrete 
equivalent of a Rayleigh-Ritz method. The Rayleigh-Ritz coefficients for the first four 
eigenfunctions with respect to the first 10 constrained basis functions are shown in 
Table 



5. Conclusions. The successful solution od a series of initial-, boundary- and 
inner-value problems, with excellent results, demonstrates the general applicability of 
the proposed matrix framework to the solution of ODEs. The generic formulation pof 
the solution method as a least squares problem is very powerful, since it enables the 
application of the methods to many classes of problems including inverse problems. 

In the case of initial value problems the least squares solution ensures that the 
solution has no implicit direction of solution, i.e., errors do not accumulate as the 
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$(^)l 


*(:c)2 


'i>(^)3 


^{x)^ 


Co 


-0.99640 


-0.06818 


-0.04144 


-0.00376 


C\ 


0.08434 


-0.93235 


-0.33425 


-0.07821 


C2 


0.00763 


0.34853 


-0.85504 


-0.36674 


C3 


-0.00317 


-0.06619 


0.39012 


-0.82201 


C4 


-0.00041 


-0.01070 


0.04490 


-0.41659 


C5 


0.00010 


0.00334 


-0.02068 


0.04126 


C6 


-0.00024 


-0.00767 


0.02352 


-0.08609 


C7 


0.00016 


0.00522 


-0.01475 


0.02897 


C8 


-0.00005 


-0.00172 


0.00487 


-0.00615 


cg 


0.00002 


0.00053 


-0.00157 


0.00340 



Table 4.2 

The discrete Raleigh-Ritz coefficients for the first four eigenf unctions with respect to the first 10 
constrained basis functions Be for the cantilever with additional simple support (see Figure (^.Tpj. 



computation proceeds. Also the application of the framework to a selection of Sturni- 
Liouville problems has delivered results comparable with those delivered by dedicated 
Sturm-Liouville solvers. 

The key issues in this paper which led to this success are: 

1. A Lanczos process with complete reorthogonalization is used to synthesize the 
polynomial basis functions. This ensures highly accurate polynomial basis for 
the computation. 

2. A correct definition of the local differentiating matrix with consistent degree 
of approximation over the complete support. This ensures the possibility of 
correctly estimating differentials at the boundary: essential for boundary and 
initial value problems. 

3. The formulation of a method of generating orthonormal homogeneous admis- 
sible functions from constraints. The matrix containing these basis functions 
is ortho-normal, yielding optimal behavior in terms of error propagation. 
This enables the implementation of a discrete equivalent of the Rayleigh-Ritz 
method. 

4. The formulation of the solution of the ODE as a least squares approximation; 
this ensure that there is no accumulation of errors. 
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A MATRIX FRAMEWORK FOR THE SOLUTION OF ODEs: 
INITIAL-, BOUNDARY-, AND INNER- VALUE PROBLEMS 

MATTHEW HARKERt AND PAUL O'LEARYt 

Abstract. A matrix framework is presented for the solution of ODEs, including initial-, bound- 
ary and inner-value problems. The framework enables the solution of the ODEs for arbitrary nodes. 
CO , There are four key issues involved in the formulation of the framework: the use of a Lanczos pro- 

cess with complete reorthogonalization for the synthesis of discrete orthonormal polynomials (DOP) 
C " 3 i orthogonal over arbitrary nodes within the unit circle on the complex plane; a consistent definition 

^\1 ' of a local differentiating matrix which implements a uniform degree of approximation over the com- 

plete support — this is particularly important for initial and boundary value problems; a method 
of computing a set of constraints as a constraining matrix and a method to generate orthonormal 
^^ , admissible functions from the constraints and a DOP matrix; the formulation of the solution to the 

'^l^ 1 ODEs as a least squares problem. The computation of the solution is a direct matrix method. The 

worst case maximum number of computations required to obtain the solution is known a-priori. This 
makes the method, by definition, suitable for real-time applications. 

The functionality of the framework is demonstrated using a selection of initial value problems, 
Sturm-Liouville problems and a classical Engineering boundary value problem. The framework is, 
however, generally formulated and is applicable to countless differential equation problems. 
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Key words. ODEs, Boundary value problems, initial value problems, inner value problems, 
Sturm Liouville, discrete orthogonal polynomials, differentiating matrix. 
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1. Introduction. There are a number of papers in which the Taylor Matrix 
is used to compute solutions to differential equations [?, ?]. These methods use 
►^ ' the known analytical relationship between the coefficients s of a Taylor polynomial 

^SJ ' E^nd those of its derivatives s to compute a differentiating matrix D. The matrix D 

together with the matrix of basis functions arranged as the columns of the matrix B 
en , are used to compute numerical solutions to the differential equations. The method of 

the Taylor matrix was also extended to the computation of fractional derivatives [?] . 



■"nI" i The problem associated with this approach is that the computation of the numerical 



solutions requires the inversion of the Vandermonde matrix, a process which is known 
to be numerically unstable, and dependent on the degree and node placement. The 
advantage of the Taylor approach lies in its ability to yield a solution for arbitrary 
nodes. 

^^ I A Chebyshev matrix approach was presented by Sezer [?] . The approach is fun- 

^ . damentally the same as for the Taylor matrix, whereby the Chebyshev polynomials 

are used as an alternative to the geometric polynomials. The main restriction associ- 
ated with the Chebvshev polynomial approach is that the numerical solution to the 
differential equations is restricted to the locations of the Chebyshev points; this lacks 
the generality needed for many differential equations and applications. 

Podlubny introduced a matrix approach to discrete fractional calculus [?] and 
later extended this work to partial fractional differential calculus [?, ?]. Triangular 
strip matrices play a central role in the work; they are used to perform integration. 
They implement the integration from a lower to an upper bound (or vice versa), 
whereby the errors accumulate as the integration proceeds. This poses a problem if 
inverse problems are addressed, since it gives the solution an implicit direction and a 
different accumulation of errors if the problem is solved from lower to upper bound 
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or from upper to lower. Furthermore, it is assumed that the initial value is zero. This 
makes the method unsuitable for arbitrary boundary conditions. An early source of 
this formulation was proposed by Courant et al. [?], (a later English translation of 
the paper is available [?]). 

A matrix solution specific to Sturm-Liouville problems was presented by Amodio [?] 
The method is specifically restricted to Sturm-Liouville problems; furthermore, it only 
supports solutions on regularly spaced nodes. The results are correspondingly modest 
for problems where the Chebyshev points yield better solutions, e.g., in the solution 
of the truncated hydrogen equation. A number of matrix approaches based on the 
Numerov method, and modifications of this method, have also been presented [?] for 
the solution of Sturm-Liouville problems, however, these methods can not be extended 
to ODEs in general. 

In this paper we formulate a general matrix framework for the solution of ordinary 
differential equations, with arbitrary initial-, boundary-, or inner values, the main 
contributions of the paper are: 

1. The proposal of a consistent framework of matrices and solution approaches 
which can be applied to initial-, boundary-, and inner-value problems; 

2. The implementation of new approached to the synthesis of discrete orthonor- 
mal basis functions, with and without weighting; 

3. Generating differentiating matrices which are of constant degree of approx- 
imation over the complete support. It is particularly important that the 
degree of approximation is consistent at the ends of the support if initial and 
boundary value problems are to be solved satisfactorily; 

4. The derivation of a means of synthesizing constrained basis functions which 
form orthonormal matrices. This basis functions span the space of all solu- 
tions which fulfil the constraints. They can be used as admissible functions 
in a discrete equivalent of a Rayleigh-Ritz method; 

5. The formulation of the solution of the ODEs as least squares approximations. 
In this manner there is no accumulation of errors. 

This paper is structured as follows: In Section ^ the framework for the gener- 
ation of all the matrices required to formulate differential equations as matrix linear 
differential operators is presented. Section ([3]) presents the approach to discretization 
of the differential equations and their solution as a least squares minimization is pre- 
sented. The required conditions for a unique solution are derived and two solution 
approaches are presented: a direct solution in the case of a unique solution and the 
implementation of a discrete Rayleigh-Ritz method for eigenvalue/eigenvector solu- 
tions, e.g., as encountered in the solution of Sturm-Liouville problems.. Finally, in 
Section Q the performance of the proposed framework is tested with a series of initial- 
value problems, Sturm-Liouville problems and a classical Engineering boundary value 
problem. 

2. Algebraic Framework. In this section we derive the structure and methods 
for the synthesis of all matrices required for the discretization and solution of ordinary 
differential equations. 

2.1. Quality Measure for Basis Functions. An objective measure for the 
quality of a set of basis functions is required if the sources of numerical error are 
to be determined and the best synthesis method is to be selected. In this paper 
continuous polynomials are considered which form orthogonal bases when evaluated 
over a discrete measure. The basis functions 6^, i.e., the polynomials evaluated at 
discrete points, can be concatenated to form a matrix, B = [bi . . .6„]. The discrete 



orthogonal polynomials (DOP) are characterized by the relationship, 

B'^WB = I, (2.1) 

where W is the weighting matrix. The Gram matrix is defined as G = B""" W B. 
Consequently, the orthogonal complement G-*" = I — B""" W B = should be a matrix 
containing only zeros. However, this is not the case, due to the loss of orthogonality 
in the three term relationship resulting from numerical errors. These numerical errors 
determine the quality of the basis functions and for which we require a measure. The 
determinant of G has in the past been used as a measure for the quality eg — det G 
of the basis functions. However, this measure does not yield stable estimates [?, 
Chapter 2, Sec. 2.7.3]. We propose the Frobenius norm of G-*^ as an error measure, 
i.e., €p = IIG^IIj;^, this is the sum of the square of all errors w.r.t. the orthogonality 
of the basis functions, ep > 0. This is a posteriori measure, i.e., we compute the 
basis functions and then determine their quality. Wilkinson [?] points out that a- 
priory prediction of error bounds yield unreliable results and a posteriori analysis is 
preferred. The numerical results obtained for different synthesis procedures can be 
found in Section (14.11). 



2.2. Numerically Stable Synthesis of Basis Functions and their Deriva- 
tives. Gram [?] proposed what is now known as the Gram-Schmidt orthogonalization 
process to generate polynomials [?]. The Gram-Schmidt process is, however, numer- 
ically unstable [?, Chapter 5] and errors accumulate as the number of integrations 
increases, i.e., with increasing polynomial degree. This precludes the synthesis of 
polynomials of higher degree with this method. Considerable research has been per- 
formed on discrete polynomials and their synthesis [?, ?, ?, ?, ?, ?, ?, ?]. The research 
was primarily in conjunction with the computation of moments for image processing. 
None of these papers present a method which is capable of synthesizing discrete or- 
thogonal polynomials of high quality for arbitrary nodes located within the unit circle 
on the complex plane. 

Here it is proposed to synthesize the polynomial basis functions using a Lanczos 
process with complete reorthogonalization [?, Chapter 9, p. 482],[?]. The procedure 
can be summarized as follows: Given a vector x of n nodes with mean x, i.e., the 
points at which the differential equation is to be solved: first compute the two basis 
functions bo, bi and initialize the matrix of basis functions B, 

bo = l/V^. bi= ,1 '^~!',| and B = [bo,6i]. (2.2) 

The remaining polynomials are synthesized by repeatedly performing the follow- 
ing computations: 

1. Compute the polynomial of the next higher degreqj, 

6„=6io6„_i; (2.3) 

2. perform a complete reorthogonalization, 

b„ = 5„ - B B^ b„ (2.4) 

= {I-BBT}6„ (2.5) 



^The symbol o represents the Hadamard product. 



by projection onto the orthogonal complement of all previously synthesized 
polynomials. It is important to note that the reorthogonalization is w.r.t. to 
the complete set of basis functions, not just the previous polynomial. 

3. Normalize the vector, 

bn = jr^, (2.6) 

4. and augment the matrix of basis functions, 

B=[B,b„]. (2.7) 

This procedure yields a set of orthonormal polynomials from a set of arbitrary nodes 
located within the unit circle on the complex plane. Although in [?] the Lanczos 
process is used to compute discrete orthogonal polynomials, the authors seem to have 
overseen the possibility (necessity) of using complete reorthogonalization at each step 
of the polynomial synthesis. 

By taking the derivative of the recurrence relationship w.r.t. x, we obtain the 
equations required to simultaneously synthesize the differentials of the polynomials. 
This procedure appears in [?] for the Legcndre and Chebyshev polynomials. Here the 
method is generalized to the synthesis of polynomials from arbitrary nodes. With 
this, the synthesis procedure delivers a set of orthonormal basis functions B and their 
derivatives B. 

2.3. Weighted Basis Functions. A set of discrete basis functions in matrix 
form, Bu, are orthogonal with respect to a weighting matrix W if, 

b;i;wb^ = I. (2.8) 

In the case of a weighting function w{x) the weighting matrix is given by W = 
diag {w {x i) .. .w{xn)}- Given a set of orthonormal basis functions B and a posi- 
tive definite weighting matrix W, there exists a set of weighted basis functions B„, 
such that B„ = B U, whereby U is a full rank upper triangular matrix. Substituting 
into Equation (|2.8p yields, 

U^B^WBU^I. (2.9) 

Since U is full rank, we may invert it to obtain, 

B'^WB^ U^'^U"^ (2.10) 

The Cholesky decomposition chol {A} of a matrix exists and is unique such that 
A = G G""" if A is real positive definite. The matrix G is a full rank lower triangular 
matrix. Consequently, the Cholesky decomposition chollB""" WB} exists if W is real 
positive definite, since B is orthonormal. Applying the decomposition yields, 

B^WB = GG'^ = U"'^U"^ (2.11) 

The sought matrix U is clearly given by, 

U = G"^. (2.12) 

With this the weighted basis functions are fully defined. The condition number of the 
basis functions depends soley on the condition number of the weighting matrix W. In 
the case where the weighting matrix is derived from a weighting function w{x)^ the 
condition number is determined by the extreme values of w{x). 



2.4. Differentiating Matrices. There are both global [?, ?, ?, ?, ?] and lo- 
cal [?, ?, ?, ?, ?] approaches to computing discrete estimates for derivatives. Global 
methods proposed in the past have used the known relationship between the coef- 
ficients of a polynomial and the coefficients of the derivative of the polynomial to 
compute a differentiating matrix. 

The computation of a differentiating matrix from polynomial bases proceeds as 
follows: The spectrum of the signal y with respect to the basis functions B is computed 
as, 



B+y. 



(2.13) 



For example, B may be the Vandcrmonde matrix; this is case with Taylor methods [?, 
?, ?]. The relationship between the spectrum s and the spectrum of the derivatives 
is given by. 



whereby 
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Consequently, 



y=BMB+y = Dy. 



(2.15) 



That is, the differentiating matrix is computed as, D = B M B+. In the case of 
the Taylor (Vandcrmonde) matrix this involves computing the pseudo-inverse of the 
Vandermonde matrix: with all the associated numerical problems. In the case of the 
Chebyshev polynomials [?, ?] B+ = B""" and a different matrix M is required, see [?] 
for details. The method is not appropriate if arbitrary nodes are required, e.g. this 
may be required if the framework is to be used to solve the problems associated with 
monitoring mechanical structures [?] . The advantage of Global methods is that they 
deliver a differentiating matrix which is valid for the complete support. 

The solution chosen here is to compute D form the basis functions and their 
derivatives, i.e., given B and B an appropriate derivative operator, D, should have the 
property that. 



Post-multiplying by B"^ yields 



DB = B. 



Db = DBB^ ^ BB' 



(2.16) 



(2.17) 



If the basis function set is complete, i.e., B B""" then the above equation yields the 
differentiating matrix directly, 



D = BB'r. 



(2.18) 



This computation is valid for arbitrary nodes. If a truncated, i.e., an incomplete, set 
of basis functions is used then B B'^ is the projection onto the the basis functions B 
and Db is then a regularizing differentiating matrix. The matrix D^ can be applied 
to the computation of estimates for derivatives in the presence of noise. 
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Fig. 2.1. Rank deficiency of the differentiating matrix D as a function of the degree of a Gram 
polynomial, when D is synthesized using Equation i2.18[ l. 



A differentiating matrix should be rank-1 deficient; it should have the constant 
vector as its null space, i.e., Dla = 0. The properties of the D are, however, depen- 
dent on the nodes being used, e.g. the Chebyshev nodes permit global differentiating 
matrices for very high degrees. With other sets of nodes the condition number of a 
differentiating matrix can increase with the degree of the polynomial being used. At 
some point the matrix starts to have additional null spaces, which are associated with 
numerical errors occuring due to insufficient numerical precision, this effect is shown 
in Figure (|2.ip for the Gram polynomials. 

Once the condition number of a global differentiating matrix has degenerated 
below an acceptable level, it becomes necessary to compute local approximations [?, ?]. 
Courant [?, ?] proposed, in 1927, using both forward and backward differences to 
compute estimates for the first derivative. The method has also been used in [?, ?] 
to this end. More commonly the tri-diagonal matrix, shown here for 6 points, 
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is used to compute a discrete local estimate for the first differentiao. This operator 
is only of degree d — 2 accurate in the core of the approximation and at both ends 
of the support only of degree d = I accurate. This makes this discrete operator 
unsuitable for the computation of derivatives at the end of the support, as is required 
for BVPs and IVPs. It is also not suitable for systems whose solutions are locally of 
degree higher than d> 2. Furthermore, this operators assumes equally spaced nodes. 
In [?] higher order finite difference schemes are proposed with end-point formulas. 
However, only equidistant spaced nodes are considered. For example, the appropriate 
three point operator for the Gram Dq^^ and for the Ghebyshev Dc,3 nodes, are. 
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^This is the matrix embedded in the Matlab function gradient. 



and 
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both computed for n — 6 points in the interval — 1 < x < 1. Note the three points 
formulas at both ends of the support. 

In keeping with the formulation of the basis functions for arbitrary nodes: the 
method for local differential approximation is also formulated here for arbitrary nodes. 
A generalized formulation of local differentiating matrix requires the vector a; of n 
arbitrarily placed nodes, the support length Is and the degree d of the approximation. 
Only odd support lengths are considered here, to avoid the need for forward and 
backward formulas. It is convenient to define the support length Is = 2ws + 1 in 
terms of the half-width Ws- The vector x of nodes is segmented into m — n ~ 2ws 
overlapping segments, for each segment. 



s{i) = x{i — Ws : i + Ws) Vi G [ws + l,n — Ws 



(2.22) 



a local set of basis functions B^ and derivatives of the basis functions Bs are com- 
puted. Then the differentiating matrix associated with the segment is determined 
Dg = Bs BJ. The first and last segments yield the end-point formulas as required. 
The remaining segment yields the required central formula of coefficients to locally 
approximate the derivative. Clearly, for the inner-segments it is only necessary to 
compute the center row vector of the local differentiating operator D^. The use of 
approximating or interpolating polynomials leads to the generation of differentiating 
matrices with and without regularization respectively. The Wilkinson diagram for the 
general structure of a local differentiating matrix D^ is shown in Equation (j2.23|) for 
the example oi Is = 5 and n = 10. The specific entries in the matrix are a function of 
the spacing of the nodes. 
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(2.23) 



All computations of the local derivative are of length Is and of constant approximation 
degree da — 2ws over the complete support. This is important if derivatives are to be 
computed at the ends of the support; furthermore, errors at the end of the support 
associated with inconsistent approximations will propagate through the entire solution 
when D is being used in the solution of differential equations. This procedure proposed 
here delivers a local differentiating matrix for arbitrary nodes. 



2.5. Defining Constraints. In Section (|2.4I) it was shown that a discrete ap- 
proximation to differentiation can be computed as a linear matrix operator. Con- 
sequently, both differential and integral constraints are linear. In the framework 
proposed here, a constraint is implemented by restricting a linear combination c^ y 
of the solution vector y to have a scalar value d, i.e., 

c^y = d. (2.24) 

This is a very general mechanism, since any constraining function can be implemented 
at a point Xi for which a linear n point expansion around this point exists. To give an 
example, consider the C^ continuous periodicity constraint y{0) = y(l), y(0) = y{l) 
and y{Q) = 2/(1): given the differentiating matrix D and D^, the three constraints can 
be formulated as: 

[l,0,...,0,-l]?/ = cTy = 0, (2.25) 

{D(l,:)-D(end,:)}y = cjy = 0, (2.26) 

{D2(1, :) - D^iend, :)} y = cj y = 0. (2.27) 

Given a set of m constraints, the constraining vectors c^ are concatenated to form 
the matrix C = [ci . . . c„i] and the corresponding scalars di form the vector d = 
[di . . .dm], so that, 

Cy = d. (2.28) 

2.6. Homogeneously Constrained Admissible Functions. Starting from 
a set of basis functions B such that B'^ W B — I , we wish to derive a method of 
synthesizing a set of constrained basis functions Be which fulfil the conditions: 

BJWBc = I, C'^Bc = and B^ = BX, (2.29) 

i.e., the constrained basis functions form an orthonormal basis set with respect to 
the weighting matrix W. If B is orthonormal, i.e., B""" B = I then so is Be. The 
constrained basis functions fulfil the homogeneous constraints defined by C. If B is 
complete then it spans the complete nx n space, given p — rank (C), i.e., the number 
of independent constraints, B^ is of dimension n x (n — p) and spans the complete 
space in which the constraints are fulfilled. Consequently, all possible vectors y which 
fulfil the constraints are given by, 

y=Bca (2.30) 

where a is an n — p vector. 

A solution to the task of determining X was presented in [?]; however, a more 
succinct derivation is provided here. The conditions from Equation (|2.29p require, 

C'^BX = (2.31) 

and with this X must lie in the null space of C''" B. Applying QR decomposition to 
BT C yields, 

QR = B'^C, (2.32) 

and consequently, 

X^QR = (2.33) 



The matrices Q and R are partitioned according to the span and null space of B""" C, 



Q = [Q„Q„] and R 



(2.34) 



with Rs of dimension p x p. The n x p matrix Q^ forms a basis set for the span and 
the n X [n — p) matrix Q„ forms a basis set for the null space of B C. Consequently, 

X^Q, =0 and (X^ Qn)'^ WX^ Q„ = I. (2.35) 

Now applying an RQ decomposition to Q„ yields, 

RQ„ = Q„. (2.36) 

R is orthonormal, since both Q„ and Q„ are by definition orthonormal. Now, selecting 
X = R yields X""" R Q„ = Q„, and with this all the conditions from Equation (j2.29p are 
fulfilled. The matrix X being orthonormal ensures that Be fulfils the same orthonormal 
condition as does B. Furthermore, X has an implicit partitioning. 



X = 



Xi 

X2 



(2.37) 



whereby, Xi is a p x (n — p) block matrix and X2 is a (n — p) x (n — p) upper triangular 
matrix. This structure ensures that the number of roots in the constrained basis 
functions Be is ordered in the same manner as in B. 

3. Discretizing and Solving Ordinary Differential Equations. In the pre- 
vious section all the matrices required for the discretization of ordinary differential 
equations were derived. In this section the discretization of initial-, boundary- and 
inner value problems is presented together with the associated methods of solving the 
resulting matrix equations. 

3.1. Initial Value Problems. In this paper we are considering the solution of 
linear ordinary differential equations with constant or variable coefficients, they can 
in general be formulated as, 

pu{x)y^''\x) . . . + pi{x) y^^\x) + po{x) y{x) = g{x) (3.1) 

to which a set of k constraints are required to ensure a unique solution. The term 
yW(^x) represents the k^^ derivative of y{x). Given the matrices derived previously, 
the discretization of Equation (|3.ip is direct and simple, each term pk{x)y'^^\x) is 
dicreteized as follows: The matrix Pfc is formed such that Pfc = diag {pk{x)}, whereby 
Pk{x) is the vector of values obtained by evaluating the function Pn{x) at the vector of 
points x\ the term y'^^\x) is discretized as D*^ y, i.e., the fc**^ power of D, which is the 
differentiating matrix derived in Section (|2.4p . Summarizing, each term is discretized 
as follows, 

Pk{x)y'^''\x)-^PkD''y. (3.2) 

and the vector g — g{x). Applying this to all terms in Equation p.ip yields, 

PfcD*^^y... + PiDy + Poy = g (3.3) 
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The matrix equivalent ol the linear differential operator L is now defined as, 

L ^ Pfe D^ . . + Pi D + Po, (3.4) 

and the set of k constraints are implemented as defined in Section (12.51) , yielding 

L y = g given C^ y = d. (3.5) 

the matrix C has the dimension n x k. 

A unique solution to the ODE exists only if 



rank 



L 



= n (3.6) 



i.e., the linear differential operator and the constraints must form a full rank system 
of equations. There are many Engineering application where this is not the case, e.g. 
the equations for the vibration of a beam, and Sturm-Liouville problems. A different 
solution approach is proposed for this class of problems in Section p.2[) . 



3.1.1. Solution as a constrained least squares problem. The formulation 
of determining y from Equation (|3.5|) as the solution of a least squares minimization 
problem yields, 

min||Ly-g||2 given C y = d. (3.7) 

y 

This is the well known problem of least squares with equality constraints (LSE). 
Efficient and accurate solutions can be found in [?, Chapter 12]. This method will 
yield solutions for ODEs with consistent constraints and a least squares solution in the 
case of over-constrained systems and perturbed systems. It is not a suitable approach 
for Sturm-Liouville type problems. 

The worst case number of floating point operations (FLOPS) required to perform 
the computation is known a-priori. This, by definition, makes the method suitable 
for real time applications. 

3.1.2. Spectral Reqularization. Spectral regularization is introduced here to 
limit the number of zeros in the basis functions and with this to reduce the errors 
associated with aliasing. Assuming y can be sufficiently accurately approximated by 
a series of r orthonormal basis functions, we may write, 

y = B,a, (3.8) 

whereby B^ — B(:, 1 . . . r). Now defining L,. = L B,. and C,. = B^ C, and substituting 
into Equation (|3.7p yields, 

min ||L,. a — g||2 given Cj a = d, (3.9) 

a 

whereby the series coefficients a are to be determined. In addition to introducing 
regularization, the truncated basis functions also reduce the size of the LS problem 
to be solved. 
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3.1.3. Solution of Homogeneously Constrained IVPs. Homogeneously Con- 
strained IVPs for a special subclass of problems for which there is a particularly simple 
solution. Let the solution y be a linear combination of a set of constrained basis func- 
tions, i.e., y = BcOt, which fulfil the homogeneous constraints C""" Be = associated 
with the IVP. Equation (13.71) now simplifies to the unconstrained least squares prob- 
lem, 



min||LBca-g||2. (3.10) 



The solution of which is. 



a = {LBJ+ff (3.11) 

since null {L Be} = if a unique solution exists. Consequently, 

y = Be{LBJ+g (3.12) 

3.2. Sturm-Liouville and Boundary Value Problems. A Sturm-Liouville 
problem is a second order ODE with the following structure. 



_d_ 
dx 






+ g{x)y = Xw{x)y, (3.13) 



in the finite interval xi < x < a;„, where p{x), g{x) and w{x) are real-valued strictly 
positive. Additionally there are two boundary conditions which are most commonly 
formulated as, 

aiy{xi)+a2v{xi) ^Q, (3-14) 

hiv{x2) + b2y{x2)=Q. (3.15) 

There are some important properties of Sturm-Liouville equations [?] which must be 
considered when implementing a discrete solution: 

1. All eigenvalues are real and there is no largest eigenvalue, i.e., there are an 
infinite number of eigenvalues and A™ — >■ cx) as to — >■ oo. Given a set of n 
discrete points x there can theoretically only be n eigenvalues; 

2. The m}^ eigenf unction has ra zeros on the interval a < a; < 6. However, 
given n points (samples) only functions with a maximum of n/2 zeros can be 
discribed without aliasing. Consider the Sturm-Liouville equation y — Xy = 
with the constraints y{0) — and y{Tr) — 0. This equation is known to have 
the eigenfunctions ^„i (x) = v2 sin (to, tt x) . Consequently, a discrete solution 
can only model the first n/2 eigenpairs correctly. 

3. The eigenfunctions are orthogonal with respect to the weighting function 
w{x), i.e., /^ w{x)<i>i{x)<i>j{x) = S{i,j). ^^ 

The general Sturm-Liouville problem formulated in Equation (I3.13|) with its cor- 
responding boundary conditions can be discritized directly as, 

{DPD-G} y = -AWy given C"^ y = 0. (3.16) 

whereby, P = diag{p(a;)}, G = dia.g{g{x)} and W — dia.g{w{x)}. A direct solution 
of this equation will, however, yield unstable results due to aliasing. 

We now introduce a set of weighted and constrained basis functions B^, which fulfil 
the orthogonality condition B^ W Buj = I and boundary conditions C"^ B = 0. These 
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basis functions are admissible functions for the Sturm-Liouville problem. The number 
of zeros in the basis functions increases from left to right in the matrix. The number 
of zeros in the admissible functions is limited, so as to avoid aliasing, by truncating 
to the first k = n/2 basis functions, i.e., B^ = Bi„(:,l : k). The eigenfunctions 
are now found as linear combinations of these admissible functions, i.e., y = B^ ct. 
Substituting this into Equation p.l6|) yields, 

{DPD-G} BaQ; = -AWBaa. (3.17) 

Pre-multiplying both sides by Bj now yields, 

Bj {DPD-G} B„a = -Aa. (3.18) 

since Bj\N Ba = I. Now defining La — Bj {D P D + G} Ba yields a standard eigen- 
vector problem, 

{la +X\} a = 0, (3.19) 

y=Baa. (3.20) 

Solving Equation (j3.19p for the eigenvalues A^ and the eigenvectors cti, then back 
substituting a^ into Equation p.20p yields the desired eigenfunctions. 

It is important and interesting to note the the matrix L^ is of dimension n/2 x 
n/2, in contrast to the original matrix L = D P D — G which is of dimension n x n. 
Consequently, dealing with the aliasing has also reduced the size of the eigenvalue 
problem to be solved. In the worst case an eigen-decomposition is of complexitj|j 
between 0(n^) and 0{n^). The improvement in speed is then in the range of a factor 
of 4 to 8, while simultaneously improving the accuracy of the solution. However, 
some of the computation gains are spent on additional pre- and post-calculations. 
A consequence of Equation (j3.20|l is that the matrix of eigenvectors a, contains the 
spectrum of the eigenfunctions with respect to the basis functions used, i.e., the 
Rayleigh-Ritz coefficients. 

4. Performance Testing. In this section a selection of examples are presented 
to demonstrate the functionality of the proposed methodQ. 

4.1. Quality of Basis Functions. The first test addresses the quality of basis 
functions, since these form the basis for all subsequent calculations. The following 
polynomials are compared: a set of Gram polynomials generated using Gram-Schmidt 
orthogonalization [?] ; a set of Chebyshev polynomials generated using the recurrence 
relationship [?] ; a Vandermonde matrix and a set of polynomials synthesized using 
the method proposed in this paper. The Frobenius norm of the projection onto the 
orthogonal complement of the Gram matrix is used as an estimate of the total error. 
The number of significant digits is then estimated to be d = — logiQ{ep). Two com- 



putations were performed: Figure (4.1(a) I shows the result for complete polynomial 



sets, i.e., the degree d = n — 1 where n is the number of nodes; Figure (4.1(b) ) is for 
a fixed number of nodes n — 1000 and the degree of the polynomial is progressively 
increased. The results shown in Figure (|4.ip indicate that the algorithm presented in 



•^Indeed there are more efficient algorithms; however, their complexity depends on the structure 
of the matrix and the distance between the eigenvalues. Consequently, no general statements can be 
made about these methods. 

^A MATLAB toolbox DOPbox is available at , http://www.mathworks.de/matlabcentral/fileexchange/41250l 
which implements all the methods proposed in this paper; furthermore, the 
code for a ll the following exampl es in also provided in the toolbox: ODEbox 

|http://www. mathworks.de/matlabcentral/fileexchange/ 1 
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Cliebysliev 



Gram 



Vandermonde 



(a) Complete polynomial basis sets, i.e., d = n- 
1 where n is the number of nodes. 



Chebytihev 




200 300 

Degree 



(b) Number of nodes n = 1000, the degree of 
the polynomial is increasing. 



Fig. 4.1. Estimate of the number of significant digits for different polynomials as a function of 
degree. DOP refers to the discrete orthonormal polynomial synthesis as proposed in this paper. 



this paper generates the most stable sets of polynomials. For this reason the algorithm 
is used for the generation of all bases required in this paper. 

4.2. Initial Value Problems. In each of the following examples the analytical 
solution is compared with the numerical results commuted using the newly proposed 
method and a Runga-Kutta solution with variable step sizqf|. The ode 4 5 is used 
for comparison in all the following initial value problems. 

4.2.1. IVP Example 1. The first example is a second order initial value prob- 
lem with constant coefficients, the equation is [?], 



2/-6y-92/ = 0, with, y(0) = 10 and 2/(0) = -75. 
in the interval < x < 3 The analytical solution to this equations is, 

y{x) = 10e~3^-45a;e"3^. 



(4.1) 



(4.2) 



The analytical solution and the results of the numerical solutions are shown in Fig- 



ure (4.2(a)). Non-uniformly spaced nodes were used, with a higher density of nodes 



where y has a higher first derivative. This demonstrates the possibility of generating 
basis functions from arbitrary nodes. The residual errors, see Figure (4.2(b)), with 
the new method is 7 orders of magnitude smaller than with the ODE4 5 method. 

4.2.2. IVP Example 2. The second example is a second order differential equa- 
tion with variable coefficients, the equation is [?], 



2x^ y — xy — 2y = 0, with, y[l) — 5 and y{l) ~ 0. 
in the interval 1 < a; < 10. The analytical solution to this equations is. 



y{x) =x +-^. 
Jx 



(4.3) 



(4.4) 



The comparison of the analytical solution with the numerical solutions using the new 
method and a Runge-Kutta procedure is shown in Figure ()4.3p . The residual error 
with the new method is 3 orders of magnitude smaller that with the Runga-Kutta 
procedure. This example has demonstrated the ability of the proposed method to 
solve initial value problems with variable coefficients. 



'See the MATLAB documentation for details of the ODE45 solver used for this computation. 
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Analytical 

o New 

V Runge-Kutta 




10 

a 5 



-5 



(a) Comparison of the analytical solution, the 
numerical solutions using the new method (with 
support length Is = 13, the nodes are placed 
X = 5z^ and z has n = 85 evenly spaced points 
in the interval < z < 1) and using a Runge- 
Kutta procedure. The nodes are shown below 
the curve. 




(b) Residual errors; (top) between the analyti- 
cal solution and the new numerical procedure, 
(bottom) between the analytical solution and 
the numerical solution using the Runga-Kutta 
procedure. 



Fig. 4.2. Computation results for the IVP Example 1 given in Equation 




X 10' 




(a) Comparison of the analytical result, the nu- 
merical solutions using the new method (with 
support length l^ = 13, there are n = 73 evenly 
spaced nodes in the interval 1 < x < 10) and 
using a Runge-Kutta procedure. 



(b) Residual errors; (top) between the analyti- 
cal solution and the new numerical procedure, 
(bottom) between the analytical solution and 
the numerical solution using the Runga-Kutta 
procedure. 



Fig. 4.3. Computation results for the IVP example 2 given in Equation {4-3\ 



4.2.3. IVP Example 3. The third example [?] is a third order non-homogeneous 
differential equation with constant coefficients, the equation is, 

y +3i/ + 3?) + y ^aOe^"^ with, y(0) = 3, 2;(0) = -3, and i/(0) = -47 

(4.5) 
in the interval < x <8. The analytical solution to this equation is, 



y{x) = (3-25a;2 + 5x3)e" 



(4.6) 



The comparison of the analytical solution with the numerical solutions using the 
new method and a Runge-Kutta procedure is shown in Figure (|4.4p . Once again 
the residual error with the new method is orders of magnitude smaller that with the 
Runga-Kutta procedure. 
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V Runge-Kutta 
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(a) Comparison of the analytical solution with, 
the numerical solutions using the new method 
(with support length Is = 13, there are n = 73 
evenly spaced nodes in the interval < a:: < 8) 
and a Runge-Kutta procedure. 




(b) Residual errors; (top) between the analyti- 
cal solution and the new numerical procedure, 
(bottom) between the analytical solution and 
the numerical solution using the Runga-Kutta 
procedure. 



Fig. 4.4. Computation results for the IVP Example 2 given in Equation J4-5|). 



4.3. Sturm-Liouville Problems. The test package for Sturm-Liouville Solvers [?] 
has been used as a source of test cases in this sectior^. 

4.3.1. Sturm-Liouville Example 1. The simplest Sturm-Liouville problem is 
chosen as the first example, since the analytical solution is known. This enables the 
investigation of the stability of the numerical computation, i.e., how many eigenvalues 
can be computed to a given accuracy. It is the equation of a vibrating string, 



y + Xy = 0, given y{0) — 0, and y(7r) = 



(4.7) 



in the interval < x < tt. The analytical solution yields the analytical eigenvalues Xk 
and eigenfunctions yk, 



Xk — k, yk ~ sin kx for fc = 1 , 



(4.8) 



The discrete solution has been computed on the interval < a; < tt sampled at the 
corresponding Chebyshev points; however, scaled so that the first and last points lie 
exactly at and tt respectively. For a given set of n points m = n/2 constrained 
basis functions are computed which fulfil the boundary values. These are the admis- 
sible functions used in what is essentially a discrete equivalent of the Rayleigh-Ritz 
method. Two different computations ni = 100 and n2 = 1000, have been performed 
to investigate the behavior of the solution with respect to the number of points used. 
A support length Is = 13 was used during the generation of the differentiating ma- 



trix. The results can be seen in Figure (4.5(a) I and (4.5(b) I respectively. The method 
returns the coefficients of the series of admissible functions required to generate each 
eigenfunction. The matrix of these coefficients is denoted by R. We use the matrix 
S — logj^Q {abs(R)} as a visual representation for the spectrum in the figures, since at 
S(i, j) ~ —16 the numerical resolution of computation environment is reached. This 
makes a visual recognition of when the spectrum fails simple. 

With Til ~ 100 approximately ki — 28 and for 712 = 1000 approximately ^2 = 280 
eigenvalues could be computed with a relative error smaller than 0.1%. This result 



''At this point we feel it is important to note that the framework presented here is generally 
applicable to the solution of ODEs in general and is not a dedicated Sturm-Liuville solver. The use 
of Sturm-Liouville problems as a test cases is to demonstrate this generality. 
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is significantly better than all previously reported results with matrix methods for 
Sturm-liouville problems [?] . This confirms the numerical stability of the approach. 
It also verifies that the number of eigenfunctions which can be computed to a given 
accuracy scales linearly with the number of nodes used. 




10 20 30 40 

Eigenvector number 71 

(a) Solution for n = 100. 




100 200 300 400 500 

Eigenvector ntunber n 

(b) Solution for n = 1000. 



Fig. 4.5. Solution of the Sturm- Liouville problem corresponding to the vibrating string for 
points. Top: spectrum of the eigenfunctions with respect to the admissible functions {logio{S). 
Bottom: The relative error in the eigenvalue Aj. in %. 

4.3.2. Sturm-Liouville Example 2. The second example is a Mathieu differ- 
ential equation; we have taken this example from [?, Problem 2, with r = 25]. This 
equation arises in the vibration of elliptical membranes. We have chosen this prob- 
lem because it is known to produce a pair of closely located eigenvalues; this should 
enable the test of the resolution of the eigenvalues computed using the proposed 
method. Secondly, the solution to the equation has no known analytical form, this 
makes numerical solutions particularly valuable. The Mathieu differential equation 
is, 



y -I- 2 r cos (2x) y = A y, with y(0) = 0, and y(n) = 



(4.9) 



in the interval < x < tt. A Chebyshev distribution of n = 1000 nodes are used 
covering the complete interval, a support length Ig = 13 and m ~ 500 basis functions 
were used for the computation. The result of the numerical computation can be seen 



in Figures (4.6(a) I and (4.6(b)). The pair of expected eigenvalues are computed as, 
Ai = -2.131489^ + 01 and A2 = -2.131486£;-t- 01. 

The results of these computations are slightly more difficult to interpret abso- 
lutely since the analytical results are not given. It can be said that the expected 
pair of eigenvalues have been found. Secondly, from the nature of the Rayleigh-Ritz 



spectrum and the corresponding computed eigenvalues, see Figure (4.6(b)), suggest 
that approximately the first k — 350 eigenvalue are correctly computed. 

4.3.3. Sturm-Liouville Example 3. This is the truncated hydrogen equation, 
taken from [?, Problem 4]. This is an example of s singular Sturm-Liouville equation 
with a limit point non-oscillatory (LPN) end point. Although only one constraint is 
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100 200 300 400 500 
Eigenvector number n 

(a) The first two eigenfunctions of tfie Mathieu (b) Top: Rayleigh-Ritz Spectrum of the Math- 
differential equation, for r = —25; note the neg- icu differential equation. Bottom: the eigenval- 
ative sign for the value of r. ues. 

Fig. 4.6. Solution of the Mathieu differential equation. It is important to note that the coeffi- 
cient r is negative. 



available the equation is well conditioned, since g{x) — >■ (X) as x — >■ 0. Highly accurate 
results for some eigenvalue are available for this equation [?]. These values are used 
to evaluate the accuracy of the computation method presented here. The differential 
equation is, 



y 



1 



y = Xy, with y{0) = LPN, and y(lOOO) = (4.10) 



in the interval < x < 1000. 

The computation was performed using n — 1000 Chebyshev points on the range 
< X < 1000; further, m = 500 basis functions were used, whereby only one constraint 
is applied, i.e., at x = 1000 and a support length of Is = 13. The result of the 
computation are presented in Table (|4.1|) . The known eigenvalues, i.e., Ao, Ag, A17 
and A18 are comparable up the the 10*'\ S***, 6*^ and 5**^ significant digits respectively. 
This indicates a high degree of accuracy, particularly considering that this is a general 
framework for differential equations and not a dedicated Sturm-Liouvalle solver. In [?] 
difficulties with oscillations at the right end point were observed for eigenvalues A^ 
when fc > 8, these difficulties are not observed with the methods proposed here. 





new method 


known value [?] 


Rcl. Error 


Ao- 


-6.2499999978£;-02 


-6.2500000000£; - 02 


3.4874503285£;-10 


Ag- 


-2.0661156136£;-03 


-2.0661157025£;-03 


4.3009091823£;-08 


Al7 = 


-2.5757218232£;-04 


-2.5757359232£;-04 


5.4741446402£;-06 


A18 — 


2.8740937561£;-05 


2.8739013100^;- 05 


-6.6963370220S-05 



Table 4.1 

Table of eigenvalues for the truncated hydrogen equation, 
obtained using the new method with the known results [7] 



Comparing the numerical results 
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4.4. Boundary Value Problem Example 1. The last test is a classical Engi- 
neering boundary value problem. A cantilever with additional simple support forcing 
the constraint 2/(0.8) = is shown in Figure (I4.7P : this is an example of an inner 
constraint. Furthermore, the system is over constrained, since there are 5 constraints 
placed on a 4*^ order differential equation. It demonstrates the ability of the pro- 
posed framework to solve problems with arbitrarily placed constraints and to solve 
over-constrained systems. The first two admissible and eigenf unctions are shown in 
Figures (4.8(a)) and (4.8(b) I respectively. This computation was performed with 



1/(0) = 



r(i) = o 



1/(0) = 



y(0.8)=0 A i/(l)=0 



Fig. 4.7. Cantilever with an additional support representing an inner-value problem. 
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(a) The first three admissible functions for the 
cantilever with and additional simple support. 



(b) The first three eigenfunctions. 



Fig. 4.8. Admissible functions and eigenfunctions for the cantilever with additional simple 
support shown in Figure j4- 



n — 1001 points evenly distributed along the interval < a; < 1, a total of r = 500 
basis functions were used and a support length Is — 13 was used for the computa- 
tion of the local derivatives. The method presented for this problem is a discrete 
equivalent of a Rayleigh-Ritz method. The Rayleigh-Ritz coefficients for the first four 
eigenfunctions with respect to the first 10 constrained basis functions are shown in 
Table 



5. Conclusions. The successful solution od a series of initial-, boundary- and 
inner-value problems, with excellent results, demonstrates the general applicability of 
the proposed matrix framework to the solution of ODEs. The generic formulation pof 
the solution method as a least squares problem is very powerful, since it enables the 
application of the methods to many classes of problems including inverse problems. 

In the case of initial value problems the least squares solution ensures that the 
solution has no implicit direction of solution, i.e., errors do not accumulate as the 
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$(^)l 


*(:c)2 


'i>(^)3 


^{x)^ 


Co 


-0.99640 


-0.06818 


-0.04144 


-0.00376 


C\ 


0.08434 


-0.93235 


-0.33425 


-0.07821 


C2 


0.00763 


0.34853 


-0.85504 


-0.36674 


C3 


-0.00317 


-0.06619 


0.39012 


-0.82201 


C4 


-0.00041 


-0.01070 


0.04490 


-0.41659 


C5 


0.00010 


0.00334 


-0.02068 


0.04126 


C6 


-0.00024 


-0.00767 


0.02352 


-0.08609 


C7 


0.00016 


0.00522 


-0.01475 


0.02897 


C8 


-0.00005 


-0.00172 


0.00487 


-0.00615 


cg 


0.00002 


0.00053 


-0.00157 


0.00340 



Table 4.2 

The discrete Raleigh-Ritz coefficients for the first four eigenf unctions with respect to the first 10 
constrained basis functions Be for the cantilever with additional simple support (see Figure (^.Tpj. 



computation proceeds. Also the application of the framework to a selection of Sturm- 
Liouville problems has delivered results comparable with those delivered by dedicated 
Sturm-Liouville solvers. 

The key issues in this paper which led to this success are: 

1. A Lanczos process with complete reorthogonalization is used to synthesize the 
polynomial basis functions. This ensures highly accurate polynomial basis for 
the computation. 

2. A correct definition of the local differentiating matrix with consistent degree 
of approximation over the complete support. This ensures the possibility of 
correctly estimating differentials at the boundary: essential for boundary and 
initial value problems. 

3. The formulation of a method of generating orthonormal homogeneous admis- 
sible functions from constraints. The matrix containing these basis functions 
is ortho-normal, yielding optimal behavior in terms of error propagation. 
This enables the implementation of a discrete equivalent of the Rayleigh-Ritz 
method. 

4. The formulation of the solution of the ODE as a least squares approximation; 
this ensure that there is no accumulation of errors. 



